2021-11-17

Aufgabe 1

Frage: Gibt es einen Gradienten entlang des Breiten- und Längengrads?

Ausgangsgrafik in DS1

Testwahl und Hypothese

Testwahl

  • Gerichtete Beziehung: Breiten- bzw. Längengrad beeinflussen den Fang (cpue_log), aber nicht anders rum
  • Es geht erstmal primär um eine lineare Beziehung (was anderes habt Ihr bisher nicht gelernt): lineares Regressionsmodell

Hypothese:

  1. H1 bei Breitengrad: Es gibt einen (bzw. keinen bei H0) linearen Effekt des Breitengrads auf den logarithmierten CPUE, sprich der (logarithmierte) CPUE ändert sich linear mit einer Änderung des Breitengrads.
  2. H1 bei Längengrad: Es gibt einen (bzw. keinen bei H0) linearen Effekt des Längengrads auf den logarithmierten CPUE, sprich der (logarithmierte) CPUE ändert sich linear mit einer Änderung des Längengrads.

Annahmen der linearen Regression

Vorab (a priori) testbar:

  • Normalverteilung der Antwortvariablen cpue_log

Anschließend (posteriori) zu prüfen (bei Modellvalidierung):

  • Normalverteilung der Residuen
  • Homogenität der Residuen
  • Keine einflussreichen Datenpunkte (Cooks Distanz < 1)

Normalität von cpue_log

Grafische Überprüfung

  • H0: Daten sind normal verteilt.
  • H1: Daten sind nicht normal verteilt.

Kolmogorov-Smirnoff-Test

# Da evtl. zeitliche und räumliche 
# Abhängigkeit gegeben ist, ist 
# dieser Test geeigneter:
y <- dat$cpue_log
ks.test(y, "pnorm", 
  mean = mean(y), sd = sd(y))
    One-sample Kolmogorov-Smirnov test

data:  y
D = 0.052766, p-value = 0.6398
alternative hypothesis: two-sided
  • \(\Rightarrow\) Die Annahme der Normalverteilung ist beim logarithmierten CPUE gegeben (Kolmogorov-Smirnoff-Test, D = 0.053, p = 0.64).

Lineare Regression

Beispiel Längengrad

mod <- lm(cpue_log ~ long, dat)

\(\Rightarrow\) Vor der Betrachtung des numerischen Outputs kommt aber die Modellvalidierung!

Modellvalidierung

Grafisch

Modellvalidierung

Normalitätstest der Residuen

y <- residuals(mod)
ks.test(y, "pnorm", mean = mean(y), sd = sd(y))
    One-sample Kolmogorov-Smirnov test

data:  y
D = 0.058788, p-value = 0.5006
alternative hypothesis: two-sided
  • \(\Rightarrow\) Die Annahmen scheinen erfüllt zu sein. Es kann nun der numerische Output betrachtet werden.

Numerischer Output

summary(mod)
Call:
lm(formula = cpue_log ~ long, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4206 -0.7396 -0.0511  0.9667  3.0094 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.48228    0.52837   17.95   <2e-16 ***
long        -0.38599    0.03582  -10.78   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.262 on 196 degrees of freedom
Multiple R-squared:  0.3721,    Adjusted R-squared:  0.3689 
F-statistic: 116.1 on 1 and 196 DF,  p-value: < 2.2e-16

Der Output in den rote Boxen ist relevant!

Zusammenfassung

Die H1 kann bestätigt werden. Die lineare Regressionsanalyse zeigt, dass es einen signifikanten, negativen Einfluss des Längengrads auf den logarithmierten CPUE der Scholle gibt (p < 0.001). Mit jedem Grad Richtung Osten, nimmt der logarithmierte CPUE um 0.386 Einheiten ab:

\[ln(CPUE) = 9.48 -0.386 * Längengrad\]

Ingesamt lässt sich 37% der CPUE Variabilität über den Längengrad erklären.

Aufgabe 2

Frage: Unterscheidet sich der Fang zwischen den ICES subdivisions?

Datenexploration

dat %>%
  group_by(area) %>%
  count()
# A tibble: 7 × 2
# Groups:   area [7]
   area     n
  <int> <int>
1    21    17
2    22    29
3    23     3
4    24    50
5    25    81
6    26    16
7    28     2

Entfernung von SD 23 und 28

dats <- dat %>%
  filter(area %in% c(21,22,24,25,26))
dats$area <- factor(dats$area)

Frage: Unterscheidet sich der Fang zwischen den ICES subdivisions?

Ausgangsgrafik in DS1

Testwahl und Hypothese

Testwahl

  • Gruppenunterschiede und > 2 Gruppen: 1-faktorielle ANOVA oder Kruskal-Wallis Test

Hypothese

  • H0: Es gibt keinen Unterschied im durchschnittlichen (ln-transformierten) CPUE zwischen den fünf ICES subdivisions.
  • H1: Es gibt einen Unterschied im durchschnittlichen (ln-transformierten) CPUE zwischen den fünf ICES subdivisions.

Annahmen der linearen Regression

Vorab testbar:

  • Normalverteilung der Antwortvariablen cpue_log IN JEDER GRUPPE!
  • Varianzhomogenität von cpue_log ZWISCHEN DEN GRUPPEN!

Posteriori zu prüfen (bei Modellvalidierung):

  • Normalverteilung der Residuen ÜBER ALLE GRUPPEN ZUSAMMEN
  • Homogenität der Residuen ZWISCHEN DEN GRUPPEN!
  • Keine einflussreichen Datenpunkte (Cooks Distanz < 1)

Normalität von cpue_log

SD21

Da evtl. eine Abhängigkeit zwischen den Daten existiert, ist der Kolmogorov-Smirnoff-Test besser. Bei kleinem N (wie hier) wäre aber Shapiro-Wilk Test besser:

Beispiel SD21 (Kattegat) 

y <- dats$cpue_log[dats$area == "21"]
ks.test(y, "pnorm", mean = mean(y), 
  sd = sd(y))
    One-sample Kolmogorov-Smirnov test

data:  y
D = 0.13831, p-value = 0.901
alternative hypothesis: two-sided
shapiro.test(y)
    Shapiro-Wilk normality test

data:  y
W = 0.95358, p-value = 0.5156
  • \(\Rightarrow\) Beide Tests deuten auf eine Normalverteilung der Stichprobe aus SD21 hin (D = 0.14, p = 0.9 bzw. W = 0.96, p = 0.5).

Normalität von cpue_log

Alle SDs

Shapiro-Wilk

dats %>% 
group_by(area) %>%
  summarise(p_vals = shapiro.test(cpue_log)$p.value)
# A tibble: 5 × 2
  area   p_vals
  <fct>   <dbl>
1 21    0.516  
2 22    0.712  
3 24    0.183  
4 25    0.835  
5 26    0.00493

Bei SD26 muss die H0 abgelehnt werden; hier sind die Daten nicht normal verteilt.

Varianzhomogenität von cpue_log

  • Grundsätzlich passt der Bartlett- oder Levene-Test.
  • Wenn aber eine der Gruppen nicht normal verteilt ist (z.B. SD26), ist der Levene-Test geeigneter:
car::leveneTest(y = dats$cpue_log, group = dats$area)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   4  1.4042 0.2342
      188
  • \(\Rightarrow\) Varianzen sind homogen. Da die Annahmen der Homogenität und Normalität (außer bei SD26, aber das ist vernachlässignar, wenn Daten homogen sind) erfüllt sind, kann die ANOVA durchgeführt werden.

ANOVA

mod <- aov(cpue_log ~ area, dats)

\(\Rightarrow\) Vor der Betrachtung des numerischen Outputs kommt aber die Modellvalidierung!

Modellvalidierung

Grafisch

Numerischer Output

  • Es gibt keine zu einflussreichen Werte und die Residuen sind auch normal verteilt und homogen.
  • Es kann daher der numerische Output betrachtet werden:
summary(mod)
             Df Sum Sq Mean Sq F value   Pr(>F)    
area          4  160.7   40.18   24.32 3.14e-16 ***
Residuals   188  310.6    1.65                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post hoc Test für paarweisen Vergleich

TukeyHSD(mod)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = cpue_log ~ area, data = dats)

$area
              diff       lwr        upr     p adj
22-21 -0.007549236 -1.089087  1.0739889 1.0000000
24-21 -0.519785941 -1.513850  0.4742780 0.6024862
25-21 -1.371264612 -2.315831 -0.4266983 0.0008614
26-21 -3.404010057 -4.637282 -2.1707380 0.0000000
24-22 -0.512236705 -1.338685  0.3142120 0.4320102
25-22 -1.363715376 -2.129914 -0.5975170 0.0000200
26-22 -3.396460821 -4.499101 -2.2938211 0.0000000
25-24 -0.851478671 -1.488266 -0.2146912 0.0027570
26-24 -2.884224116 -3.901206 -1.8672419 0.0000000
26-25 -2.032745445 -3.001402 -1.0640891 0.0000003

Zusammenfassung

Die H0 kann abgelehnt werden, es gibt einen signifikanten Unterschied zwischen den Gruppen (1-faktorielle ANOVA: \(F_{4, 188} = 24.32\), p < 0.0001).

Der Tukey’s HSD Test zeigt, dass sich hauptsächlich SD25 und 26 von den anderen 3 signifikant unterscheiden. Hier ist der logarithmierte CPUE im Durchschnitt geringer, insbesondere in SD26 (s. Abb. links).

Aufgabe 3

Frage: Was könnten mögliche Gründe für ein Muster sein? Könnte die Tiefe eine Rolle spielen, oder die Temperatur oder der Salzgehalt?

Ausgangsgrafik in DS1

Testwahl und Hypothese

Testwahl

  • Gerichtete Beziehung: Umweltparameter beeinflussen den Fang (cpue_log), aber nicht anders rum.
  • Es geht erstmal primär um lineare Beziehung (was anderes habt Ihr nicht gelernt): lineares Regressionsmodell

Hypothese:

  1. H1 bei Tiefe: Es gibt einen (bzw. keinen bei H0) linearen Effekt der Tiefe auf den logarithmierten CPUE, sprich der logarithmierte CPUE ändert sich linear mit einer Änderung der Tiefe.
  2. …

Annahmen der linearen Regression

siehe Aufgabe 1

Vorab (a priori) testbar:

  • Normalverteilung der Antwortvariablen cpue_log

Anschließend (posteriori) zu prüfen (bei Modellvalidierung):

  • Normalverteilung der Residuen
  • Homogenität der Residuen
  • Keine einflussreichen Datenpunkte (Cooks Distanz < 1)

Multiple Regression - Kollinearität zwischen Xs

Zusätzliche Annahme VORAB zu testen:

VIF < 3? 

mod <- lm(cpue_log ~ depth + temp + sal, dat)
car::vif(mod)
   depth     temp      sal 
1.591708 1.257096 1.576777
  • \(\Rightarrow\) Alle unabhängigen Variablen haben einen VIF < 3, also gibt es keine Kollinearität!

Multiples Regressionsmodell

mod <- lm(cpue_log ~ depth + temp + sal, dat)
  • \(\Rightarrow\) Vor der Modellauswahl und Betrachtung des numerischen Outputs kommt aber die Modellvalidierung!

Modellvalidierung

Grafisch

Automatische Modellauswahl basierend auf AIC

step()

step(mod)
Start:  AIC=123.17
cpue_log ~ depth + temp + sal

        Df Sum of Sq    RSS    AIC
- depth  1     0.036 354.26 121.19
<none>               354.22 123.17
- temp   1    32.204 386.42 138.40
- sal    1    50.280 404.50 147.45

Step:  AIC=121.19
cpue_log ~ temp + sal

       Df Sum of Sq    RSS    AIC
<none>              354.26 121.19
- temp  1    39.639 393.90 140.19
- sal   1    73.227 427.48 156.39
Call:
lm(formula = cpue_log ~ temp + sal, data = dat)

Coefficients:
(Intercept)         temp          sal  
   -4.16317      1.68793      0.08493

Finales Modell nur mit temp und sal

fin_mod <- lm(cpue_log ~  temp + sal, dat)

Nochmal Modellvalidierung

Numerischer Output

summary(fin_mod)
Call:
lm(formula = cpue_log ~ temp + sal, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6267 -0.8119  0.0496  0.9677  4.0506 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.16317    1.40369  -2.966   0.0034 ** 
temp         1.68793    0.36135   4.671 5.57e-06 ***
sal          0.08493    0.01338   6.349 1.49e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.348 on 195 degrees of freedom
Multiple R-squared:  0.2875,    Adjusted R-squared:  0.2802 
F-statistic: 39.34 on 2 and 195 DF,  p-value: 4.433e-15

Der Output in den rote Boxen ist relevant!

Zusammenfassung

Text

\[ln(CPUE) = -4.16 + 1.688 * Temperatur + 0.085 * Salinität\]

Das finale Regressionsmodell zeigt, dass der logarithmierte CPUE linear mit einer Zunahme der Temperatur und Salinität steigt (jeweils p < 0.0001). Allerdings kann das Modell nur 28% der Gesamtvariabilität des Fangs erklären, ein großer Teil der räumlichem Variabiltät muss daher durch andere Umweltbedingungen zu erklären sein.

Grafik

  • Wenn mehr als eine Kovariate im linearen Regressionsmodell sind, kann nicht mehr einfach mit der predict(model) Funktion die Regressionslinie eingezeichnet werden.
  • Auch die geom_smooth() Funktion ist nicht geeignet, da sich diese immer auf eine einfache Regression bezieht!

Grafische Zusammenfassung

Temperatur

library(modelr)
df_temp <- data_grid(data = dat,
  temp = seq_range(temp, 20),
  sal = mean(sal) # keep constant
)
pred_temp <- predict(fin_mod, newdata = df_temp, 
  se.fit = TRUE)
df_temp$pred <- pred_temp$fit
df_temp$ci_low <- pred_temp$fit - pred_temp$se.fit
df_temp$ci_up <- pred_temp$fit + pred_temp$se.fit

ggplot(dat, aes(x = temp)) +
  geom_point(aes(y = cpue_log)) +
  # add regression line
  geom_line(data = df_temp, aes(y = pred), 
    col = "blue") +
  # add confidence intervals
  geom_ribbon(data = df_temp, aes(
    ymin = ci_low, ymax = ci_up), 
    fill = "grey50", alpha = 0.2)

Grafische Zusammenfassung

Salinität

df_sal <- data_grid(data = dat,
  temp = mean(temp), # keep constant
  sal = seq_range(sal, 20) 
)
pred_sal <- predict(fin_mod, newdata = df_sal, 
  se.fit = TRUE)
df_sal$pred <- pred_sal$fit
df_sal$ci_low <- pred_sal$fit - pred_sal$se.fit
df_sal$ci_up <- pred_sal$fit + pred_sal$se.fit

ggplot(dat, aes(x = sal)) +
  geom_point(aes(y = cpue_log)) +
  # add regression line
  geom_line(data = df_sal, aes(y = pred), 
    col = "blue") +
  # add confidence intervals
  geom_ribbon(data = df_sal, aes(
    ymin = ci_low, ymax = ci_up), 
    fill = "grey50", alpha = 0.2)